Predicción del número de vehículos registrados en el sistema de tránsito nacional

1. Introducción

En Colombia, el Registro Único Nacional de Tránsito (RUNT), funciona como una gran base de datos centralizada en la que se almacena información sobre todos los vehículos del país, los conductores, los seguros de estos, las infracciones entre muchos más. El presente reporte asume el ejercicio de crear un modelo estadístico para predecir el número de vehículos registrados diariamente en el RUNT durante el año 2018.

2. Datos

La base de datos inicial, registro_autos_entrenamiento, en formato .XLSX, contiene los registros diarios de autos en el RUNT, desde las fechas 1/01/2012 hasta 31/12/2017. Se usaron los software Excel y R para el procesamiento de la base de datos. Inicialmente, la base de datos contiene 2192 observaciones y 2 variables. Las variables Fecha y Unidades. Una vista previa de las primeras observaciones de la base de datos es la siguiente:

datos <- read_excel("registros_autos_entrenamiento.xlsx") #cargar datos

head(datos, n=5)

Con el objetivo de plantear un modelo más interesante, se propone la creación de cinco nuevas variables que permitirán desarrollar un modelo predictivo que arroje más información en sus resultados teniendo en cuenta el año, mes, días del mes, días de la semana y si el día tuvo o no alguna festividad. Las variables creadas son:

YEAR: año del registro. DIA: día del mes del registro.
MES: mes del registro. DIA_SEMANA: día del registro (lunes a domingo) FESTIVIDAD: indica si el registro se realizó en alguna de las siguientes fechas especiales: a_nuevo: 1 de enero de cada año. día_habil: día en donde no es festivo no hay festividades. festivo: día festivo. sem_santa: día de semana santa, de sábado a domingo. madres: día de las madres. brujas: día de las brujas, 31 de octubre. navidad: 25 de diciembre. a_viejo: 31 de diciembre.

Todas las nuevas variables exceptuando FESTIVIDAD se obtienen de la variable Fecha. En el caso de la variable FESTIVIDAD, se obtiene a través una base de datos externa que se adiciona a la base de análisis original y abarca los días feriados en Colombia desde 2012 hasta 2018.

Para efectos practicos en la sección de modelo predictivo, se exportan dos bases de datos:

La primera, con 7 variables y 2192 observaciones que van desde el año 2012 hasta el año 2017.

La segunda, con 6 variables y 2575 observaciones que van desde el año 2012 hasta el año 2018. Esta base de datos no contiene la variable UNIDADES debido a que se borró intencionalmente para luego rellenarla con la predicción del modelo.

3. Análisis descriptivo

plot_ly(data = descriptivo, x = ~FECHA, y = ~UNIDADES,
        type = "scatter", mode = "lines", split = ~YEAR,
        line = list(width = 1)) %>%
  layout(title = 'Registros de vehiculos diarios en el RUNT de 2012-2017',
        xaxis = list(title = "Año"),
        yaxis = list(title = "Número de vehiculos"))

Aprecia que gradualmente van disminuyendo los registros diarios, se registran cada vez más vehículos los últimos días de cada mes durante todo el año, el promedio diario de vehículos registrados parece mantenerse ya que el último día del mes compensa la falta de registros del resto de días pero en realidad el registro de vehículos si disminuye año por año ya que esta compensación no es suficiente Además, en esta grafica se puede apreciar una gran cantidad de días en los cuales no se registran vehículos, en los cuales resaltar cada 1 de enero y 25 de diciembre.

descriptivo$DIA_SEMANA <- ordered(descriptivo$DIA_SEMANA, levels = c("lunes","martes","miercoles","jueves","viernes","sabado","domingo"))


totals5 <- descriptivo %>%
  group_by(DIA_SEMANA) %>%
  mutate(total=sum(UNIDADES))
totals5

ggplot(descriptivo, aes(fill = DIA_SEMANA, x =  DIA_SEMANA, y = UNIDADES)) +
  geom_bar(stat = "identity") +
  geom_text(data = totals5, aes( y = total, label = total), vjust = -0.5) +
  scale_y_continuous(breaks = seq( 0,450000,50000)) +
  xlab("Días") + 
  ylab("Numero de vehiculos") + 
  ggtitle("Registro de vehiculos por día de la semana") +
  theme(legend.position = "none")

Podemos apreciar que de miércoles a viernes se hace la mayoría de los registros, con el viernes como el día de la semana con mayor cantidad de vehículos registrados, se puede observar como hay menos registros en los fines de semana en comparación a los días de la semana, ademas de que el domingo es el día con menor cantidad de registros de vehículos.

totals4 <- descriptivo %>%
  group_by(MES) %>%
  mutate(total=sum(UNIDADES))


ggplot(descriptivo, aes(fill = MES, x = as.factor(MES), y = UNIDADES)) + 
  geom_bar(stat = "identity") +
  geom_text(data = totals4, aes( y = total, label = total), vjust = -0.5) +
  scale_y_continuous(breaks = seq( 0,200000,25000)) +
  scale_fill_viridis_d( option = "G") +
  xlab("Mes") +
  ylab("Numero de vehiculos") +
  ggtitle("Registro de vehiculos por mes")+
  theme(legend.position = "none")

Diciembre es el mes con mayor cantidad de registros de vehículos y enero el mes con menor cantidad de registros, ignorando estos dos datos y solo tomando del mes febrero al mes de noviembre los registros de vehículos en el RUNT están entre 130000 y 150000, variando por unos cuantos miles de vehículos

totals3 <- descriptivo %>%
  group_by(YEAR) %>%
  mutate(total=sum(UNIDADES))


ggplot(descriptivo, aes(fill = YEAR, x = as.factor(YEAR), y = UNIDADES)) + 
  geom_bar(stat = "identity") +
  geom_text(data = totals3, aes( y = total, label = total), vjust = -0.5) +
  scale_y_continuous(breaks = seq( 0,400000,25000)) +
  scale_fill_viridis_d( option = "E") +
  xlab("Año") +
  ylab("Numero de vehiculos") +
  ggtitle("Registro de vehiculos  por año")+
  theme(legend.position = "none")

El año 2014 es el año con más vehículos registrados, ignorando este dato podría decirse que cada año se registran menos vehículos.

totals2 <- descriptivo %>%
  group_by(dia_habil) %>%
  mutate(total=sum(UNIDADES))

gfes1 <- ggplot(descriptivo, aes(fill = dia_habil, x = dia_habil, y = UNIDADES)) +
  geom_bar(stat = "identity") +
  geom_text(data = totals2, aes(x = dia_habil, y = total, label = total), vjust = -0.5) +
  scale_y_continuous(breaks = seq( 0,2250000,250000)) +
  xlab("dia habil") +
  ylab("Numero de vehiculos") +
  ggtitle("registros de dias habiles y no habiles")+
  theme(legend.position = "none")

totals <- solo_tipo_festividad %>%
  group_by(FESTIVIDAD) %>%
  mutate(total=sum(UNIDADES))


gfes2 <- ggplot(solo_tipo_festividad, aes(fill = FESTIVIDAD, x = reorder(FESTIVIDAD,+UNIDADES), y = UNIDADES)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_d( option = "D") +
  geom_text(data = totals, aes(x = FESTIVIDAD, y = total, label = total), vjust = -0.5) +
  xlab("festividad") +
  ylab("Numero de vehiculos") +
  scale_y_continuous(breaks = seq( 0,30000,2500)) +
  ggtitle("Registros en las festividades") +
  theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "none")


ggarrange(gfes1,gfes2)

A comparación de los días hábiles, la suma de vehículos registrado en los días festivos y de festividades es muy poca. Los dos unicos días que no presentan ningún registro de vehículos en el RUNT de 2012 a 2017 son año nuevo(primero de enero) y navidad(25 de diciembre).

Siendo muchos menos los días festivos comunes que los días celebración de la semana santa se registran muchos más vehículos en semana santa, puede hacerse una comparación parecida con el día de madres y haloween, siendo solo un día hay una enorme diferencia en los registros siendo el día de madres 2 registros al ser día domingo y en haloween más de 11mil al ser casi siempre un día normal en semana.

4. Modelo Predictivo

Se lleva a cabo un Modelo lineal generalizado con la familia de distribución Poisson, donde la variable respuesta (variable dependiente) sería “UNIDADES”, y las variables explicativas (variables independientes) son: FECHA, YEAR, DIA, MES, DIA_SEMANA, y FESTIVIDAD.

Donde se utilizarán los años desde 2012 hasta 2016 para la creación del modelo.

Se procede a hacer lectura de la base de datos en donde se almacenarán todas las predicciones, según el modelo anterior.

Se muestran las primeras 10 observaciones de las predicciones dadas. Almacenadas en la variable NRO_TOTAL_VEHICULOS, que representa el número de vehículos registrados.

prediccion <- predict(object = lm2, newdata = Base_prediccion,
                          type = "response")

prediccion_diaria<- Base_prediccion %>% 
  mutate(Nro_Vehiculos = round(prediccion,0))

diario_20_02 <- prediccion_diaria %>%
  group_by(FECHA, YEAR, DIA, MES, DIA_SEMANA, FESTIVIDAD) %>%
  mutate(NRO_TOTAL_VEHICULOS=Nro_Vehiculos)

diario_20_02$NRO_TOTAL_VEHICULOS <- as.integer(diario_20_02$NRO_TOTAL_VEHICULOS)

head(diario_20_02, 10)

Luego, se observa el MSE, MAE y \(R^2\) para los datos de entrenamiento (que posee los días comprendidos entre el 01/01/2012 y el 31/12/2016).

#Entrenamos con 2012-2016

datos_02 <- subset(los_datos, (YEAR != '2017'))

y_train <- round(predict(lm2, newdata= datos_02, type="response"))
y_actual <- datos_02$UNIDADES
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm7_tmse, lm7_tmae, lm7_r2)

## [1] "MSE: 55163.222770, MAE: 150.340449, R2 Score: 0.818973"

Por consiguiente, se observa el MSE, MAE y \(R^2\) para los datos de validación (que posee los días comprendidos entre el 01/01/2017 y el 31/12/2017).

#Validación con 2017

datos_vl <- subset(los_datos, (YEAR == '2017'))

y_train <- round(predict(lm2, newdata= datos_vl, type="response"))
y_actual <- datos_vl$UNIDADES
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm7_tmse, lm7_tmae, lm7_r2)

## [1] "MSE: 82417.772603, MAE: 166.572603, R2 Score: 0.704797"

También se calcula el MSE, MAE y \(R^2\) para el primer semestre de 2017.

# Primer semestre del 2017

datos_vl <- subset(los_datos, (YEAR == '2017'))
sem2017<- datos_vl[c(1:181),]

y_train <- round(predict(lm2, newdata= sem2017, type="response"))
y_actual <- sem2017$UNIDADES
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm7_tmse, lm7_tmae, lm7_r2)

## [1] "MSE: 68703.745856, MAE: 145.447514, R2 Score: 0.705406"

Luego, el \(R^2\) del conjunto de validación que comprende el primer semestre del 2018 (01/01/2018 y el 30/06/2018), no se calcula ya que para el 2018 no se obtienen observaciones reales para dicho año.

A continuación, se aprecian las primeras 10 y las últimas 10 observaciones del primer archivo plano (01/01/2012 y el 31/12/2017), que contiene 1827 observaciones con dos variables.

datos_16 <- subset(diario_20_02, (YEAR != '2018'))
dat <- datos_16[,-2]
dat01 <- dat[,-2]
dat02 <- dat01[,-2]
dat03 <- dat02[,-2]
datos_primer_archivo_plano <- dat03[,-2]
head(datos_primer_archivo_plano, 10)

tail(datos_primer_archivo_plano, 10)

write.csv(datos_primer_archivo_plano, "primer_archivo_plano.csv")

Después, se pueden observar las primeras 10 y las últimas 10 observaciones del segundo archivo plano, el cual contiene el primer semestre de 2018 (01/01/2018 y el 30/06/2018). que posee 181 observaciones con dos variables.

datos_17 <- subset(diario_20_02, (YEAR == '2018'))
sem2018 <- datos_17[c(1:181),]
dat04 <- sem2018[,-2]
dat05 <- dat04[,-2]
dat06 <- dat05[,-2]
dat07 <- dat06[,-2]
datos_segundo_archivo_plano <- dat07[,-2]
head(datos_segundo_archivo_plano, 10)

tail(datos_segundo_archivo_plano, 10)

write.csv(datos_segundo_archivo_plano, "segundo_archivo_plano.csv")

Nota: Los dos archivos planos en formato .csv se pueden observar en el repositorio de Github.